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ABSTRACT 

We present a method, based on Bayesian statistics, to fit the dust emission parameters in the far-infrared and 
submillimeter wavelengths. The method estimates the dust temperature and spectral emissivity index, plus their 
relationship, taking into account properly the statistical and systematic uncertainties. We test it on three sets of 
simulated sources detectable by the Herschel Space Observatory in the PACS and SPIRE spectral bands (70- 
500 /im), spanning over a wide range of dust temperatures. The simulated observations are a one-component 
Interstellar Medium, and two two-component sources, both warm (HII regions) and cold (cold clumps). We first 
define a procedure to identify the better model, then we recover the parameters of the model and measure their 
physical correlations by means of a Monte Carlo Markov Chain algorithm adopting multi-variate Gaussian 
priors. In this process we assess the reliability of the model recovery, and of parameters estimation. We 
conclude that the model and parameters are properly recovered only under certain circumstances, and that false 
models may be derived in some case. We applied the method to a set of 91 starless cold clumps in an inter-arm 
region of the Galactic Plane with low star formation activity, observed by Herschel in the Hi-GAL survey. Our 
results are consistent with a temperature independent spectral index. 
Subject headings: TBW 



1. INTRODUCTION 

Previous observations in the submillimeter (submm) and 
far-infrared domains have enabled the investigation of the 
properties of dust in a variety of environments and at different 
angular resolutions. The brightness of radiation emitted from 
a source in local thermal equilibrium has a continuum spec- 
trum which can be expressed as I\ = e\B\{T), where e\ is 
the wavelength dependent emissivity, and B\{T) is the Planck 
function corresponding to a temperature T. In optically thin 
regions, the emissivity equals the optical depth T\. In op- 
tically thick regions the emi ssivity tends to 1 In the case 
of op tically thin dust clouds dDesert et al.lll990t iDraine & Lil 
120071) . the simplest models assume that the emissivity de- 
pends on the wavelength as a power law, and the brightness 
spectrum is: 



by: 
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where Td is the dust temperature, eo is the emissivity at wave- 
length Ao, and j3 is the emissivity spectral index, with a fidu- 
cial value j3 = 2. 

More refined models describe the spectrum as due to a com- 
bination of components different in temperature and nature. 
In particular, observational data evidence a flattening of the 
spectral index towards long wavelengths (A > 500 fj,m) which 
is well re presented by the Multi component Dust Model de- 
scribed in Fink beiner et al.l(11999l) . The spectrum is then given 



marcella.veneziani@ipac.caltech.edu 

1 Infrared Processing and Analysis Center, California Institute of Tech- 
nology, Pasadena, CA 91125, USA 

2 Dipartimento di Fisica, Universita di Roma "La Sapienza", Rome, 
Italy 

3 Universite de Toulouse, UPS-OMP, IRAP, Toulouse, France 

4 CNRS, IRAP, Toulouse, France 



/A(ei ! e2,7i,r 2 ) = ei 



Ac 



Bx(Ti) + e 2 



Bx(T 2 ) 



(2) 

where e\ and e 2 are the emissivities at Ao, T\ and T 2 are the 
temperatures of the two dust components, and fa and fa are 
the spectral indices. Standard value for spectral indices can 
be s et to: fa = 1.67, fa = 2.7 0, i.e. the best fit values in model 
8 of iFinkbeiner et all Il999l) . 

More sophisticated models incorporate the effect of 
the di sordered internal structure of amorphous dust 
grains dMenv et al.l 120071) . In this case, the emission is 
characterized by a single temperature and a spectral index 
which changes with temperature. Moreover, some level 
of anticorrelation between spectral index and temperature 
is expected by la b oratory experiments (see for exam- 
ple ICoupeaud et al.l (1201 ll) ). The emission can be then 
characterized by Eq. ([T) with the extra relation 



f3=A 



(3) 



where Tq is a pivot temperature (we use Tq = 20 K), and A and 
a are parameters that characterize the inverse relation. 

Measurements from balloon based observatories and satel- 
lites indicates the existence of an anticorrelation between 
and the T^. Th e balloon-borne e xperiments P RONAOS 
dDupac et al.l2003l) and ARCHEOPS dDesert et al.l2008t) data 
sets evidenced an inverse relationship between T</ and j3 in a 
wide range of environments of the interstellar medium (ISM) 
and in cold sources. In PRONAOS data, variations in the spec- 
tral index were observed in the range 0.8-2.4 for dust tem- 
peratures between 1 1 and 80 K, whereas ARCHEOPS data 
showed that the inverse relationship is more pronounced with 
values of (3 up to 4 for temper atures down to 7 K. Moreover, 
recently Veneziani et Id] ( 120101) highlighted a similar trend an- 
alyzing Tj and (3 for eight high Galactic latitudes clouds, by 
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combining IRAS, DIRBE and WMAP data to BOOMERanG 
observations. The f3 values vary from 1 to 5 in the tempera- 
ture range 7-20 K, with a behavior similar to the ARCHEOPS 
results. The analysis of the Herschel Hi-Gal key program 
dMolinari et al.l2010allbl) has de monstrated this antic orrelation 
in the Interstellar Medium (ISM. lParadis et al.l2010h while the 
ATLA S found similar r esults at high latitude s (IBracco et ail 
12011b . On the contrary. IPaladini et all (12012b finds the same 
inverse trend on HII regions but it might be generated from 
spurious effects. 

In order to test the hypothesis of the existence of a func- 
tional dependency between the temperature and the spectral 
index, one has to properly account for spurious correlations in 
their estimated values. Due to the spectral shape in Eq. <£TJ, the 
effect of a high value of Td can be mimicked, when holding 
the intensity fixed, by a low value of /3, and viceversa. There 
exists a statistical degeneracy between these physical quanti- 
ties that must be considered while in vestigating if an intri n- 
sic physical correlation indeed exists. IShettv et ail (l2009allbl) 
report spurious inverse correlations due to a number of fac- 
tors including noise in flux measurements coupling with the 
Tj - (3 degeneracy and variations in temperature of overlap- 
ping sources along the sightline. 

Because of the last Herschel and Planck submillimeter and 
millimeter data on dust, the topic of the 7^-/3 spurious an- 
ticorrelation has been re cently addressed (see for example 
Uuvela & Ysardll2012allbl: iKellv et aLlHoTlh both on simula- 
tions and on real data. The originality of our paper is to treat 
statistical and systematic uncertainties in two different ways, 
both during the SED fitting and in the 7^-/3 relationship es- 
timate, in order to take into account their different statisti- 
cal properties. In this paper we present a method based on 
Bayesian statistics, to discriminate between models in Eq. (fl]i- 
(f2]i taking into account both the statistical and systematic er- 
rors of the input fluxes. We also describe how to distinguish 
between a spurious and a physical anticorrelation between 
the spectral index and the temperature, and to estimate the 
Td — f3 inverse relationship by properly evaluating the degen- 
eracy caused by the spectral shape and the flux uncertainties, 
both statistical and systematic, i.e. calibration uncertainty. 

We demonstrate the robustness of this method in the Her- 
schel PACS and SPIRE spectral coverage, 70fim < A < 
500/im, and for a temperature range 10 < T < 50K. We 
consider 3 cases: the diffuse ISM, warm sources (HII re- 
gions/YSOs), and cold sources (e.g. pre-stellar clumps). 

This paper is organized as follows: Sec. (f2) presents the al- 
gorithm, after a brief introduction on the temperature-spectral 
index degeneracy and Bayesian statistics; Sec. (O describes 
the simulated data set used in the analysis; Sec. and[5]re- 
port the results obtained by applying the method to the simu- 
lated observations and to cold clumps on the Galactic plane 
detected by the Herschel/Hi-GAL survey. Conclusions are 
drawn in Sec. ©. 

2. BAYESIAN METHOD FOR IR SEDS FITTING 

In this section, after a brief discussion of the temperature- 
spectral index degeneracy and of Bayesian statistics, we 
present our method to fit the SEDs taking into account both 
the statistical noise and the systematic uncertainties present in 
every dataset. The combination of Bayesian techniques and a 
good knowledge of the parameter distribution in the parame- 
ter space, allows to remove the spurious temperature spectral 
index anticorrelation, as well as the biases on the physical 
parameter estimates, generated both from random noise and 
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FIG. 1. — Normalized SEDs for a 15K, 40K and 80K source and spectral 
index /3 = 2. The blue vertical lines identify the Herschel bands chosen to 
perform our analysis. The SED are then well constrained at low temperatures 
(dotted line) and poorly constrained at high temperatures (solid line) where 
Herschel sampling covers only the RJ part of the spectrum. 



systematic effects. 

2.1. Bayesian treatment of statistical noise 

In order to study the effect of noise (i.e. statistical uncer- 
tainties) on one-component SED fitting, and in particular the 
Td - j3 degeneracy, we first simulate one ISM like and one HII 
region like source with a 20% gaussian random noise, and 
then fit for the physical parameters using in Eq. (Q3 . Since we 
are interested in Herschel-like observations, we sample the 
SEDs in the six bands A = [70, 100, 160,250,350,500] /im. 
The simulated source temperatures are Tism =15 K and 
Thh = 80K. The reason for chosing so different temperatures 
is shown in Fig. (HJ where the normalized SED of three 
sources at 15K, 40K and 80K and their sampling with the 
Herschel bands are reported. In the Herschel bands the SED 
is well measured at low temperature (dotted line) while for 
warmer sources (dashed and solid lines) we can only sample 
the RJ part of the spectrum without any information on the 
Wien side. The spectral index is constrained by the slope of 
the RJ portion of the grey body while the temperature is con- 
strained by the peak location. In our simulation, we fix the 
spectral index to 2. The 70/iin band is usually not included 
in ISM studies because is very sensitive to very small grains 
(VSG) emission. We include it in this simulation to show the 
dependence of the parameter recovery on the error bars and 
on the SED sampling. 

To fit the SED of the simulated sources and to study the 
parameter distribution in the parameter space, as well as 
their correlation , we use a Monte Car lo Markov Chain algo- 
rithm (MCMC, iLewis & Bridiell2002l) . Hereafter, we explain 
the main idea behind this approach. 

Following Bayesian statistics, the probability to have a set 
of parameters p given a set of data d is 

P(p\d)ocP(p)P(d\p). (4) 

Here P(-\-) are conditional probability densities with argu- 
ment on the right-hand side of the bar assumed to be true, 
P(p) is the a priori probability density, or prior, of the param- 
eters set p, P(d|p) (e.g. the likelihood function) is the prob- 
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ability density of dataset d given a set of parameters p and is 
described by the formula 



A'a 



J°(d|p) oc exp I -- ^ 



b=\ 



(5) 



where N\ is the number of bands, Ij, is the value of the cho- 
sen model with parameters p in the A;, band, either Eq. ([JJ or 
Eq.(|2]i and <t\ is the statistical error associated to the flux db- 
If Eq. ([T]i is chosen, then p = [eo, Tj,/3]; otherwise, if Eq. (O 
is chosen, then p = [ei, £2)71,72]. The a priori probability 
density of a given set of parameters is the probability distribu- 
tion that expresses the parameter uncertainties before the data 
are taken into account. When nothing is known about the pa- 
rameter distributions, usually wide flat priors are adopted, as 
they assign equal probabilities to all possibilities. If, on the 
contrary, the parameter statistical distribution is known from 
previous experience, i.e. previous measurements, then a more 
specific, informative prior might be adopted. 

The posterior probability P(p|d) is estimated using the 
MCMC algorithm. Given a set p,, with likelihood L, and 
posterior probability P(p, |d), the MCMC algorithm peforms a 
random walk through the allowd parameter space, generating 
a new independent set p, + i with likelihood L, + i and posterior 
probability P(p i+1 |d). This second set is accepted according 
to a rule which also guarantees a good sampling of the prob- 
ability density in a reasonable computational time. We make 
use of the Metropolis-Hastings algorithm, which is a class of 
Markov Chains in which the new set p,+i is always accepted 

A(^1,Q=^# = W(Pm) >1 (6) 



P(Pi\<D 



This guarantees the convergence to the maximum of the like- 
lihood function. If instead A(;+ l,z) < 1, the new set is ac- 
cepted with a probability proportional to the ratio A(;+ 1,/). 
An advantage of this method is to ensure a good sampling 
of the parameter distribution in parameter space. For an ap- 
plication of this technique on millimete r and sub-millimeter 
observations, see lVeneziani et all (1201 Ol) . 

Fig. (f2| shows the 68% and 95% posterior probabilities in 
the Td — fl parameter space obtained for the two cases of ISM 
(T/sm = 15K, P = 2) and HII region (T m/ = 80K, (3 = 2), re- 
spectively. The elongated and slant shape of these posterior 
probabilities indicate that the two parameters are degenerate 
in parameter space. In case of warm sources we have a strong 
constraint on j3 but a not precise determination of the temper- 
ature, as shown in the right panel. On the contrary, in the case 
of ISM sources (left panel), the whole SED is well sampled, 
leading to a smaller degeneracy between the temperature and 
the spectral index. 

The dependence of one of the two parameters on the other, 
i.e. their statistical degeneracy, can be easily demonstrated. 
Starting from a one-temperature modified black body: 



F(T b ,(3) = e 



\\~ P 2hc 2 1 
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In the R-J regime, we can approximate with: 



An 



(7) 



(8) 



which is a linear relation in logA-logF scale. Thus we ex 



pect some level of anticorrelation between T</ and /3, as in any 
linear relation there is some level of anticorrelation between 
slope and intercept. 

The more the fluxes are noisy and the SED not well sam- 
pled, the more the two parameters are correlated and diffi- 
cult to constrain. Moreover, when combining data of different 
sources to recover correlation between parameters, there i s the 
tendency to treat systematic errors as if they were statistical, 
which is not correct. Here we treat the systematic errors in a 
proper way, and use Bayesian statystics, in order to provide 
a correct reconstruction of the parameters distributions and of 
physical correlations between parameters. 

2.2. Bayesian treatment of systematic errors 

The sources simulated in Sec. d2.U have up to 20% of sta- 
tistical uncertainty on measured fluxes. Nonetheless, the flux 
uncertainty might come also from the calibration error and 
other sources of systematic error, which do not have the same 
properties as the statistical errors (i.e. instrumental noise, 
background fluctuations). In order to properly take into ac- 
count the systematic effects, we make use of a Monte Carlo 
procedure. In the following we will refer only to the system- 
atic error due to calibration uncertainty but the same proce- 
dure can be applied to any source of systematic error. 

We would like to emphasize that this treatment is particu- 
larly important in our specific case: the recovery of a rela- 
tionship among physical parameters of different observables, 
several sources in our case. The systematic error can be the 
same in all sources, as in the case of a wrong calibration of 
one band. If this is not considered, the result can be severely 
biased. 

In order to properly take into account in our analysis the 
calibration error we use a conditional inference technique. We 
consider a n-tuple of the calibration values k = (k\ ,k% , k^ x ). 
Each value kt, has a probability distribution P(kb), which we 
assume to be uniform in the -Akb < kb < Akb range. The 
values of Ak are reported in Tab. (|2) for each band. Given 
the observed dataset d, we have a joint distribution of the pa- 
rameters p for all possible configurations of k, P(p|d,k). The 
conditional result for a single source is obtained by marginal- 
ization of this probability over all the possible values of the 
calibration uncertainties, weighted by their distribution: 



P(p\d)- 



Ak 



P(d,k|p)P(p)P(k)dk 



Ak 



where 



P(d,k|p)ocexp(--£ 



b=\ 



db-kbhip) 



<Jb 



(9) 



(10) 



Eq. (|9]l is an exstention of Eq. (|4]i where both statistical and 
systematic uncertainties are taken into account. In order to 
apply this technique, we run a Monte Carlo simulation, scat- 
tering the values of k in 100 realizations. We tried different 
numbers of calibration realizations, but we could not find any 
difference in the final results above 100 steps. In each itera- 
tion p,, a calibration uncertainty vector k M is generated. The 
fluxes d of all sources are multiplied by the same set of k M . 
The fit of the SED in each iteration is pe rform ed with a Monte 
Carlo Markov Chain algorithm (see Sec. 12. U . 

The final physical parameters which better describe the 
SED of each source are then obtained by marginalizing the 
100 calibration dependent set of physical parameters over 
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FIG. 2. — 68% and 95% contours of the two-dimensional posterior probabilities of sources at T = 15 K, /? = 2 (left panel) and T = 80 K, /3 = 2 with a 20% 
statistical error on fluxes. The independent points of the Markov Chains are denser within the 68% contour, close to the maximum of the distribution, than close 
to the boundaries. The slant shape means that T and /3 are degenerate and the elongation is larger for the warm sources where Herschel bands do not constrain 
the peak of the SEDs. The lcr error bars are obtained by marginalizing the tri-dimensional posterior probability, as we have three parameters, over the remaining 
two dimensions. 



the calibration uncertainties as in Eq. (|9). The marginaliza- 
tion i s performed with the publicly available GetDist soft- 
ware (Lewis & Bridlell2002l) . 

2.3. Method description and uncertainties analysis 

The method to recover dust physical parameters focuses on 
our treatment of statistical and systematic uncertainties and 
on a good definition of the a priori probability density in the 
MCMC. The analysis is carried out according to the following 
steps. The details of each step are reported later in this section. 

1 . We make a first, fast MCMC run on the whole sample 
of sources to check the best model for each SED and to 
study the parameter distributon in the parameters space. 
Therefore, we fit the SED of each source with both Mi 
and M2, taking into account both statistical and system- 
atic unc ertai nties (i.e. calibration errors) as described 
in Sec. ( 12.21 ). At this stage, we assume uniform wide 
priors, reported in Tab. (Qj. 

2. We determine if Mj or M2 better fit the SED of a source, 
based on the comparison of the x 2 p robability den- 
sity functions, as reported in Sec. (12. 4k if both models 
give good probability, we assign the source the simpler 
Mi; if both models give bad probability we exclude the 
source from the analysis. An example of good and bad 
fits after step (1) is reported in Fig. (|3). 

3. For the sources identified with Mi, we estimate again 
the physical parameters taking into account both statis- 
tical and systematic uncertainties, as in step 1, but this 
time with a more precise definition of the "a priori" 
probability densities. Therefore, we assume multivari- 
ate Gaussian priors, estimated from the first run. This 
procedure is described in Sec. (12. 5\ . 

4. From the new physical values the temperature-sp ectra l 
index relationship is measured as reported in Sec. (12.61 . 

2.4. Model identification 

In real observations, especially of the Galactic Plane, there 
is often an overlap of sources along the line of sight or a multi- 



TABLE 1 

Priors during the model identification (step (1)) 
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1.7, 2.0,2.3,2.7 
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List of a priori probability densities imposed on parameter set during the 
model identification procedure. We chose wide uniform (U) ranges at values 
in order not to constrain the fit results and to avoid the code to converge to 
non physical values. The explored range of f3c values is also reported. 



component temperature source. It is crucial to be able to dis- 
tinguish among different models, and identify the correct one, 
before attempting any physical interpretation. 

As a first step in our analysis, we then fit the simulated 
dataset with both models in Eq. (Q]) and (0, i.e. with a single 
temperature and two-temperatures. In Mi, the fitted parame- 
ters are the emissivity (eo), the temperature T</ and the spec- 
tral index (3. We do not consider any model with more than 
two components because with the datasets available nowa- 
days it is unlikely to have enough data points to well fit up 
to 6 parameters. For the same reason, in the two components 
case, we cannot perform a fit over all the six parameters: Ti, 
T2,ei,£2,/3i,/32- We rather set the spectral indices to a fixed 
common value j3 c = /3\ = [3%, and estimate the two emissivities 
and the two temperatures. The values explore d for j3 c in model 
2 are 1 .7, 2.0, 2.3 and 2.7, consistently with Fin kbeiner et alj 
( 119991) . We then execute the fit for M2 four times, changing 
every time the spectral index which is assumed to be the same 
for both the components, i.e. not temperature dependent. To 
make sure not to constrain the results, we assume uniform 
wide a priori probability densities on every parameter. They 
are the same for all sources and are reported in Tab. [T] 

We can then easily find the model which better fits the SED 
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by using the \ 2 probability density function. Since Mi and M2 
have a different number of degrees of freedom, we compare 
the two analyses using the cumulative distribution function 
P of a x 2 distribution with the proper number of degrees of 
freedom. The ideal fit should have P ~ 0.5. P <C 0.5 indicates 
an overestimate of the error bars while P ^> 0.5 indicates that 
the fitting model is not a good model for the dataset. In one- 
component sources, we expect, in an ideal situation, P(M\) to 
be near to 0.5, and P{M2> to be larger. The opposite is true 
for two-component sources. In an ideal situation, we expect 
to have at least one of the P(Mr, fi c ) near 0.5, and P{M\) to be 
larger. 

An example of one and two-dimensional posterior proba- 
bilities of the physical parameters of an ISM-like source and 
of an HII region combined with dust is shown in Fig. (01. In 
the top line of this figure we show the posterior probabilities 
of a 20K source (ISM) fitted with a one temperature model 
(top left panel) and with a two-temperature model (top right 
panel). The convergence is reached with the one-temperature 
model because all the three parameters (temperature, spectral 
index and emissivity) are well constrained within the bound- 
ary provided. On the contrary, the fit with two components 
(top right) is not well constrained. The four parameters (two 
temperatures and two emissivities) do not converge within 
the acceptable ranges. In the bottom line of Fig. © the fit 
of a two component sources with a one and two-components 
model is shown. The source is an HII region of 50K com- 
bined with some surrounding dust at 17K. The spectral in- 
dices of both components are set to 2. In this case, the fit with 
a two-components model converges (bottom right) with the 
four parameters well defined within the chosen ranges, while 
the fit with a one-component model (bottom left) is poorly 
constrained within the physically acceptable boundaries. 

2.5. Parameter Estimates 

Since the parameters in Eqs. (G3 and (0 are not uniformly 
distributed and are also intrinsically correlated in the param- 
eter space, their "a priori" probability density P(p) (Eq. |9]l 
must take these properties into account. This is required to 
avoid biased estimates of the source characteristics. The in- 
trinsic parameter correlation in each source is due to the spec- 
tral shape describing the source emission (Eqs. [T]and|2j an d 
it doesn't have anything to do with the physical dependency 
of one parameter on the other, which can be measured from 
an ensemble of objects. The level of the intrinsic correlation, 
in our case, depends on the random noise present in the sin- 
gle SEDs, on the SED sampling within the considered bands 
and on possible source overlap along the line of sight. Since 
those effects might vary source by source, we do not perform 
a fit of the parameters all together, but analyze the parame- 
ter statistical distribution of each source independently of the 
others. The information about the parameter statistical distri- 
bution and intrinsic correlation is contained in the covariance 
matrix, estimated for each source by means of the first MCMC 
run (step (1)). The covariance matrix is then used to build the 
new a priori probability distribution, which is not flat any- 
more as in the first run, but contains specific and definite in- 
formation about the variables. The new set of priors P(p) is 
a multi-variate Gaussian, with the center p= [logeo,/3o,7b] 
in the approximate average values of the whole sample of 
sources. The covariances of the multi-variate Gaussian pri- 
ors are estimated, for each source, by broadening 10 times 
the parameters standard deviations keeping their correlations. 
This makes us sure not to constrain the final results with a too 
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FIG. 3. — One and two-dimensional physicalparameters posterior prob- 
abilities for an ISM-like (source Al, see Tab. l|2)) and an HII (source C2) 
region-like source estimated with a MCMC ran assuming wide uniform pri- 
ors. The ISM is well fitted by M t (top left) and not by M 2 (top right). As a 
consequence the posterior probabilities are well constrained within the range 
provided when considering a one-temperature model (top left) and not with 
a two-temperature model (top right). The opposite is true for the second one: 
Mi is not well constrained (bottom left) and converges to very low value 
of the spectral index while all the unknowns in M9 converge to finite 
results (bottom right). 



tight parameters range. In other words, while the center of 
the Gaussian prior are the same for the whole set of sources, 
the parameter covariances are estimated source by source and, 
as a consequence, the fit is performed source by source. An 
example of multi-variate Gaussian priors for a single ISM- 
like, cold clump-like and HII region-like source, is reported 
in Tab. (0]i and shown in Fig. ©. 

Taking all those aspects into account, the final parameters 
estimate on each source is then performed again by applying 
Eq. (0 and (TTOb but, this time, instead of assuming uniform 
distributed priors, we describe P(p) as a multi-variate Gaus- 
sian: 

iWaexp^-i^p-pyST^p-p))^ (11) 

where Ey is the covariance matrix of each source, with the 
indices i and j running over the entire set of paramers. Ey 
elements are obtained from the correlation matrix Cy by Ey = 
Cy 07 a j where ct, cr, are the standard deviations of the param- 
eters pi and pj, obtained from the first MCMC run, multiplied 
by a factor of 10 in order to broaden the prior, retaining the 
correlation properties. 

An example of C elements are reported in Tab. (0). S con- 
tains therefore the information about the distribution and cor- 
relation of parameters. Both p and E are estimated through 
the MCMC in step (1). The final likelihood is then 
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1 Nx 
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d-klj-pT 



N„ N„ 



exp 



i=l 7=1 



(pi-pdC \pj-pjY 



dk\2) 



where N p is the number of parameters and k is the vector of 
the band calibration uncertainties on extended emission. As 
already pointed out, with our method sources are fitted one 
by one in order to take into account different values of the 
covariance matrix, due to i.e. different random noise or dif- 
ferent sampling of the SEDs, and different systematics effects 
on different sources, i.e. overlap of other objects on the line 
of sight. 

If the sources in the sample have similar noise level and 
systematic characteristics, i.e. line-of-sight temperature vari- 
ations and calibration unce rtainties, one can a lso use a hierar- 
chical Bayesian procedure (Kelly et alj|201 lh . The hierarchi- 
cal model fits the SEDs of the whole sample of sources at the 
same time, assuming a global model for the overall distribu- 
tion of source parameters. 

2.6. Estimate of the Temperature-Spectral Index relationship 

As discussed in the Introduction, an inverse relation be- 
tween temperature and spectral index has been found in sev- 
eral previous studies. We model this relation by means of 
Eq. ®. 

The comparison between Mj and M2 and the estimate of 
the functional form of the parameter distributions, are the first 
steps to estimate the temperature-spectral index relationship. 
We select as possible source candidates of a temperature de- 
pendent spectral index, the ones with Mi preferred over M2 
and with P(M\) < 0.95 (2a). This last requirement is neces- 
sary because the worst the fit, the more spurious the T,/ - (3 
anticorrelation. 

If all the previous conditions are satisfied, we estimate the 
Td—/3 trend on the selected subsample of sources, taking into 
account the calibration uncertainties again with a Monte Carlo 
procedure on the calibration uncertainties. In each iteration 
/1, i.e. for each k^, we fit sources with Mi through Eq. © 
and Eq. (fTTT i and then estimate Eq. (O by means of a MCMC 
run assuming uniform wide priors on A M and a^: < A < 3, 
-2 < a < 2. After 100 iterations, we have 100 pairs of A M 
and On, We then marginalize over the calibration errors as in 
Eq. ®: 

/.Ak 

P(A,a\x)= / P(x;k\A,a)P(A,a)P(k)dk (13) 

J-Ak 

where x = (7d,/3), P(k) is the same as in Eq. ©, P(A,a) is 
assumed to be uniformly distributed and uncorrelated, and the 
likelihood of the p-th iteration is 

1^/W-A/' 

P(x;k\A,a) oc exp | 



T rf (k^) 



7b 



(14) 



where N s is the number of sources, /3(k M ) and 7rf(k M ) are the 
spectral index and temperature of source s, estimated with the 
calibration uncertainties k M and a((3k) is the statistical error on 



/3(k). One dimensional probabilities of A and a are obtained 
by marginalization over the other parameter. 

3. SIMULATED OBSERVATIONS 

Since we are interested in studying the SED fitting using 
the Herschel bands, as a baseline of our analysis we use the 
PACS and SPIRE wavelengths, angular resolutions and noise 
properties. We complement the Herschel spectral coverage 
simulating fluxes also in the MIPS (24/xm), IRAS (lOO^rn) 
and Planck-HFI (850/iin) bands. These bands are particu- 
larly helpful to constrain the SEDs of warm and cold single 
and multi-temperatures objects whose emission is not fully 
sampled by the PACS and SPIRE instruments. Depending 
on the model, we simulate one or two-temperatures modi- 
fied black body with a 1-er statistical uncertainty of a\ = 10% 
of the total flux in each band. This value arises from the 
combination of the experimental Herschel noise confusion 
limit (few percent of the total flux) combined with some non- 
negligi ble background fluct uation and it's the value also cho- 
sen bv lShettv et al.l (l2009al) in their analysis. The calibration 
uncertainties are estimated on diffuse emission. The fluxes 
are then scattered within the error bars twice: first, in a Gaus- 
sian uncorrelated fashion, to take into account random noise; 
second, in a correlated fashion, i.e. the same bands have the 
same scatter in all sources, to take into account the calibra- 
tion uncertainties. We assume a common angular resolution, 
i.e. 4', for IRAS 100/im, or 5' when Planck 850^m is also in- 
cluded. The considered bands and the experimental character- 
istics used to perform the simulations are reported in Tab. (0. 

We simulate observations of sources with different physi- 
cal conditions. The simplest case consists of one-component 
sources, described by Eq. dTJ, with a temperature around 20K 
and a spectral index which may or may not vary with tem- 
perature (cases Ai and A2, respectively). In the following, 
we will refer to the one-component model as model 1 (Mi). 
A more complex case is given by the overlap of more than 
one source along the line of sight. Indeed, one of the key 
questions of the SED fitting is the error made when approx- 
imating a multi-temperature observation with an isothermal 
model and if this approximation ca n lead to a spurious corre- 
lation among the fitted parameters (Shetty et al. 2009b). This 
could be the case when analyzing a source with more than 
one population of dust grains or more than one source along 
the line of sight. We have considered two different scenar- 
ios, one in which the two temperatures are very close (few 
degrees apart) and one in which they are very different (tens 
of degrees). A typical source of the first case is a cold dense 
pre-stellar clump, i.e. a cold clump with a temperature (Tc) 
around 10K surrounded by an envelope of dust with an av- 
erage temperature (T^) of around 15K (case B). A typical 
source of the second case might be an HII region combined 
with foreground/background diffuse dust (case C). Both the 
considered scenarios are described by Eq. (0 and, in the fol- 
lowing, we will refer to the two-component model as model 2 
(M 2 ). 

These last two cases are also chosen to investigate the effect 
of the band coverage. A good sampling of the SED peak and 
of the Rayleigh-Jeans (RJ) side is required, in order to prop- 
erly fit both temperature and spectral index and to identify 
the underlying physical model. For example, in HII regions, 
the Herschel wavelength coverage is not enough to separate 
the two components, while including the MIPS 24-fim the 
two components are properly detected. Similarly, pre-stellar 
clumps properties are better identified including Planck-HFI 
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TABLE 2 

Simulated observing experiments 
Band (fim) Experiment Beam size Ak ex , Ak p 



24 


MIPS (Carev et al. 2009) 


6" 


15% 


4% 


10% 


70 


PACS (Poglitsch et al. 2010) 


6" 


20% 


3% 


10% 


100 


IRAS-IRIS (Miville-Deschenes & Lagache 2005) 


4' 


13.5% 


13.5% 


10% 


160 


PACS (Poglitsch et al. 2010) 


11" 


20% 


4% 


10% 


250 


SPIRE (Griffin et al. 2010) 


18" 


15% 


7% 


10% 


350 


SPIRE (Griffin et al. 2010) 


25" 


15% 


7% 


10% 


500 


SPIRE (Griffin et al. 2010) 


37" 


15% 


7% 


10% 


850 


Planck-HFI (Planck Collaboration 201 la) 


5' 


3% 


3% 


10% 



Experimental characteristics used to generate the simulated sources. The reported calibration uncertainties (Afe» and Ak po j„,) refer to diffuse emission and to 
point sources, respectively. The statistical error <r is assumed to be 10%, which is the maximum value obtained in the analysis of Hi-Gal data. It includes both 
instrumental noise and background fluctuations. 



850/im band. 

The simulated sources are described below and more details 
are provided in Tab. (|3): 

[A] ISM environment, as in model 1: single component 
sources with a Gaussian temperature distribution cen- 
tered in 20K and ranging between 15K and 25K, ob- 
served with a 4' resolution in the spectral range 100- 
500/i m. We place ourselv es exactly in the same situa- 
tion as lParadis et al.l (120101) in order to test our ability to 
recover the real relationship between temperature and 
spectral index. For this purpose, we first set the spectral 
index to 2 (sources Ai) and then we use as input a tem- 
perat ure dependen t spectr al index, both anticorrelated, 
as in iParadis et al.l ( 120101) . (sources A2) and correlated 
positively correlated (sources A3). 



Interstellar Medi' 






[B] Cold clumps (CC), based on model 2: two-component 
sources consisting of an envelope at a temperature T^ 
and a core at a temperature 1c- The spectral in- 
dices of both components are set to 2. We investi- 
gate tw o options for the rel ative emissivities as ana- 
lyzed bv lShettv etal](l2009t)l) . We also study the depen- 
dence of the fit on the spectral coverage by measuring 
the emission in the range 100-500 fjm and then includ- 
ing also the 850/im flux which constrains the RJ side of 
the SEDs. 

[C] HII regions combined with foreground/background dif- 
fuse dust, of temperatures Thii and T«, respectively. 
This case is also based on model 2. The spectral in- 
dices are set to 2 and the relative emissivities are scaled 
according to the anal ysis performed on Hi-GAL data by 
iPaladini et al.l (120121) . As in the case of the cold clumps, 
we investigate the dependence of the physical parame- 
ters recovery on the spectral coverage. We then sample 
the SEDs in the 70-500/im range and then include also 
the 24/im flux to better constrain the temperature of 
the warmer component. In this analysis we assume the 
24/im emission to be dominated by Big Grains (BG) as 
recen tly found, for example, by SOFIA (Salgado et al. 
120121) . The BG thermalize with the Interstellar Radia- 
tion Field (ISR) and we can model their emission as a 
modified grey body. If, on the contrary, the 24/im flux is 
dominated by VSG, we need to remove this point from 
the analysis and consider only the 70 - 500/im range, 



FIG. 4. — One-dime nsion al projection of the multi-variate Gaussian prior 
P(p) described in Eq. Jilt for each set of sources. The parameters are re- 
ported in Tab. l|4}. 

because they don't have the same physical properties as 
BG. 

In order to compare the temperatures obtained fitting with 
Mi the SEDs of two-temperature sources with the input val- 
ues, we make use of a d e nsity- weighted temperature intro- 
duced by iDoty & Palottil (120021) . T col . According to their 
definition, the column temperature of a source described by 
Eq. (O is given by: 



Tcol : 



eiTi+e 2 T 2 
£i+£2 



(15) 



4. RESULTS 



In this section we present the results of the SED fitting, and 
the Tj- f3 relationship estimate after properly taking into ac- 
count the systematic errors of the datasets, for each type of 
source considered. 

The multi-variate Gaussian priors of one sample source for 
each set are reported in Tab. (2) and shown in Fig. (0). 

The results are summarized in Tab. ((5) and Tab. (|6). In 
Tab. (0, first column, we provide the input models; in the sec- 
ond, third and fourth columns, the percent fraction of sources 
assigned Mi and M2 as well as the excluded sources, respec- 
tively, highlighting the true positive detections; columns five 
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TABLE 3 
Simulated sources 

Input Central temperature Relative column Resolution Spectral Range 

sources (K) densities (A™) 



Interstellar Medium (Ai) 




Ta 


= 20 






2 


4' 


100-500 


Interstellar Medium (A2) 




Ta 


= 20 






2.2(T A /T T L3 


4* 


100-500 


Interstellar Medium (A3) 




Ta 


= 20 






2(T A /T )^ S 


4' 


100-500 


Cold Clump (B 1 ) 


T C = 


10, T £ = 


15, T co/ = 


12.5 


N C =N £ 


2 


5' 


100-500 


Cold Clump (B 2 ) 


T c = 


10, T £ = 


15, T co i = 


12.5 


N c =N £ 


2 


5' 


100-850 


Cold Clump (B 3 ) 


Tr = 


10, J E = 


15, T co/ = 


10.9 


N c = 10 N E 


2 


5' 


100-500 


Cold Clump (B 4 ) 


T c = 


10, T £ = 


15, Tcoi = 


10.9 


N c = 10 N £ 


2 


5' 


100-850 


HH+dust (Ci) 


Thii = 


17, T rfrf 


= 50,T CO , 


= 17.3 


N M = 100 N„„ 


2 


4' 


70-500 


HH+dust (C 2 ) 


Thii = 


17, 1 di 


= 50, T co , 


= 17.3 


Nrfrf = 100 N H „ 


2 


4' 


24-500 



Three types of simulated sources used to test our method. In each model we generate 1000 sources with a gaussian distribution in temperature. The distribution 
is contained on the value reported in the 2nd column, with a dispersion of 10%. In case A2 there is a Tj — anticorrelation in input. 7q is set to 20K. 



TABLE 4 

Priors during the final estimate (step (3)) 



the second peak, the marginalized final parameters obtained 
with M2 are badly constrained. 



Sources 


P 


P 


cr ¥ 


C 


P(P) 


ISM 

(A) 


log eo 



T rf (K) 


-4.2 
1.8 

23 


9.7 
2.5 
85.7 


C(loge ,/3)=0.985 

C(loge ,T rf )=-0.985 

C(/3,T,,)=-0.956 


MVG 
MVG 
MVG 


CC 
(B) 


log e 



T rf (K) 


-6.5 
1.8 
13 


6.1 
1.6 
54.8 


C(loge ,,9)=0.982 

C(loge ,T rf )=-0.986 

C(/3,T rf )=-0.951 


MVG 
MVG 
MVG 


HII+ISM 

(C) 


log e 



T rf (K) 


-5.5 
1.8 

23 


34.9 

9.0 

309.3 


C(loge ,,9)= 0.969 
C(loge ,T rf )= -0.976 
C(fiJ d )= -0.924 


MVG 
MVG 
MVG 



List of a priori probability densities imposed on the final paremeter estimates. 
The covariance matrix (£) elements are obtained by multiplying the correla- 
tion matrix elements (C) with the standard deviation of the two correspondent 
parameters. The values reported are for just one single source for each kind 
of simulated observations. The overall shape of the probability density is a 
multi-variate Gaussian (MVG) where we only take the positive values. 



to eight report for which f3 c value the best M2 fit is obtained. 
In Tab. (O we give the input and measured T</-/3 relationship 
obtained with our procedure ("Bayesian output") and with a 
least square method ("Least Square output"), for comparison, 
for the sources assigned Mi. The least square method is 
based on a x 2 minimization to recover the source parameters 
and their ler errors. Since this method do not perceive any 
treatment of the systematic uncertainties, besides adding them 
in quadrature to the statistical errors, it has been tested on 
simulations where only a 10% gaussian error bars are present, 
without including the additional calibration errors. The least 
square minimization is used in both stages, the SED fitting 
and the T-/3 relationship estimate. In this last step, the statisti- 
cal uncertainties adopted on data points are the ones obtained 
from the SED fitting step. In order to claim a physical an- 
ticorrelation with the Bayesian procedure, we require a 3 a 
detection of the a parameter. 

As a general comment we would like to highlight that the 
two-component model is not well constrained if the spectral 
coverage is poor. In order to have a clear detection of more 
than one component we need to have information both on the 
RJ and the Wien part of the SEDs. If the second component is 
very faint or the spectral coverage too small to clearly detect 



4. 1 . Interstellar Medium 

The considered cases (Ai, A2, A3) are shown in Fig. <(5j and 
reported in Tab. (|5)- 





■ (a) 




I. 




-0.2 -0.1 0.0 0.1 0.2 0.3 

AT../T, 



-0.2 -0.1 0.0 0.1 0.2 03 

AT A2 /T A 



n. 



(d) 



I 




-0.3-0.2-0.1-0.00.1 0.2 0.3 



-0.3-0.2-0.1-0.00.1 0.2 0.3 



(g) 





10 15 20 25 30 35 40 

T A i CK) 



10 15 20 25 30 35 40 

T A2 (K) 



(O 



-0.2 -0.1 0.0 0.1 0.2 0.3 

AT A3 /T A 



(f) 



-0.3-0.2-0.1-0.00.1 0.2 0.3 

amp:;, 




10 15 20 25 30 35 40 

T„(K) 



FIG. 5. — Comparison between output and input temperatures and spectral 
indices for sources A (panels from (a) to (f)). All the histograms are nor- 
malized to 1, Panels from (g) to (i): T-/3 relationships. Red circles: input 
values. Black dots: values recovered with the Bayesian methos, within 1-cr 
bars. Blue dots: values recovered with a least squares method. The light and 
dark purple contours identify the recovered correlation within 1-cr and 2-<r 
error, respectively. 



Their SEDs are much better approximated by Mi than by 
M2, being the two-components model vary badly constrained 
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TABLE 5 
Compared analysis 



Input 


M,(%) 


M 2 (%) 


Excluded (%) 


Pc= 1.7(%) 


A = 2.0(%) 


A=2.3(%) 


ft.=2.7(%) 


Interstellar Medium (Ai) 


100 




















Interstellar Medium (A2) 


98 


2 





2 











Cold Clump (Bi) 


100 




















Cold Clump (B 2 ) 


99 


1 











1 





Cold Clump (B 3 ) 


89 


11 








1 


9 


1 


Cold Clump (B 4 ) 


86 


14 











14 





HH+dust (Ci) 


32 


68 





1 


48 


19 





HH+dust (C 2 ) 


4 


96 





9 


33 


54 






Comparison of the different kind of sources approximated with a one-component (Mi) an d a two-component mod el (M 2 ). We run model 2 four times setting the 
spectral index f) c to the values 1.7, 2, 2.3 and 2.7 to explore the whole range indicated by Finkbeiner et al. 1 1999). The bold face values indicate the percentage 
corresponding to the right input models. 



TABLE 6 

Temperature-spectral index relationship from simulations 



Input Model 


Input 


Bayesian Output 


Input Recovery 


Least Square Output 


Input Recovery 




A 


a 


A(k ex t,o- 10 %) 


a(k ra t,oi %) 




A(o-io%) 


a(°"io%) 




Interstellar Medium (Ai) 


2.0 


0.0 


1.8±0.1 


0.0±0.1 


True 


1.84 ±0.03 


-0.8 ±0.1 


False 


Interstellar Medium ( A 2 ) 


2.2 


-1.3 


2.1 ±0.1 


-1.3 ±0.1 


True 


2.01 ±0.03 


-1.4±0.1 


True 


Interstellar Medium (A3) 


2.0 


0.5 


1.7 ±0.1 


0.4 ±0.1 


True 


1.83 ±0.03 


-0.5 ±0.1 


False 


Cold Clump (Ei) 


2.0 


0.0 


1.7±0.I 


-0.3 ±0.1 


True 


1.06±0.01 


-1.67 ±0.02 


False 


Cold Clump (B 2 ) 


2.0 


0.0 


1.7±0.1 


-0.3 ±0.1 


True 


1.42 ±0.04 


-0.7 ±0.1 


False 


Cold Clump (B 3 ) 


2.0 


0.0 


2.0 ±0.2 


-0.3 ±0.2 


True 


0.5±0.1 


-2.3 ±0.5 


False 


Cold Clump (B 4 ) 


2.0 


0.0 


2.2 ±0.2 


0.0 ±0.2 


True 


0.8±0.1 


-1.2±0.3 


False 


HII+dust (CO 


2.0 


0.0 


1.9±0.2 


-0.4 ±0.7 


True 


3.0 ±0.5 


-2.4 ±0.3 


False 


HH+dust (C 2 ) 


2.0 


0.0 















Output T,;— /3 relationship (columns 4-6) obtained with our method, based on Bayesian statistic with multi-variate Gaussian priors, and comparison with the input 
(column 2-3) and with the results obtained by fitting the same set of sources with a least square method (columns 7-8). Our simulations include a 10% statistical 
error and a 15% -20% systematic uncertainty depending on the band. Systematic uncertainties are treated only in the procedure based on the Bayesian method. 
The least square method do not allow a correct systematic error treatment, so results based on this procedure include only a 10% statistical uncertainty. For a 
discussion of these results see text. 



when the Monte Carlo analysis of the calibration uncertain- 
ties is included. As the plots of the relative differences show 
(panels from (a) to (f)), the temperatures are systematically 
shifted towards higher values of ~ 5 - 10% and this corre- 
ponds to a systematic spectral index shift towards lower val- 
ues of ~ 10-15% of the input on average. In order to un- 
derstand the origin of this bias, we run five new sets of sim- 
ulations of the Ai model, all with 10% statistical errors but 
with [0%, 5%, 10%, 15%, 20%] calibration uncertainties, 
assuming all the bands to have the same systematic uncer- 
tainty. Since we do not have a more precise knowledge of 
the distribution of the calibration errors, we can only assume 
a uniform flat prior in the considered interval. We run the 
whole pipeline, and estimate the bias on the final tempera- 
tures (AT/To) and spectral indices (A/3//3o). Results on the 
bias study are shown in Fig. ©. When a uniform wide prior 
is assumed (top panel), which is the case we apply on simula- 
tions and on real observations, the systematic shift increases 
with the calibration uncertainty, in a non-uniform way for the 
two parameters: the spectral index seems to be affected more 
than the temperature, likely because the calibration variations 
affect the RJ slope of the SEDs more than the peak position. 
Nonetheless, if we assume to have a better knowledge of the 
calibration uncertainty distributions, which is not the case for 
the observations considered in the present paper, and assume 



Bias from uniform calibration Bias from gaussian calibration 

o.2 ; 1 1 1 TTV — ' 0-2 ^ s ^ ^ . ^ 




-0.2 F 1 1 T 1 1 -0.2 F 1 1 1 : 

-10 10 20 30 40 -5 5 10 IS 

% Cal % Cal 

FIG. 6. — Average bias on the final temperatures and spectral indices as a 
function of the calibration uncertainties. 



gaussian priors with standard deviations [0%, 3%, 5%, 7%, 
10%], the bias level is very low for both the parameters (the 
systematic shift is centered in few percents) and is consistent 
with zero. 

A good measurement of this bias would require a bet- 
ter knowledge of the uncertainty distribution. However, our 
method takes into account this effect in the estimate of the pa- 
rameter uncertainties and, therefore, the input values and the 
recovered T— f3 are consistent with the input within la. As it 
will be clear in the following, all the parameter estimates in 
this paper are likely affected by this effect. The systematic 
shift is slightly more important when more components along 
the line of sight are assumed (cases B and C) but, again, this 
effect is included in the final error bars. 



10 



M. Veneziani et al. 



4.2. Cold clumps 

The results of the analysis on cold clumps type of sources 
is reported in Tab. (01 (cases Bl, B2, B3 and B4) and shown 
in Fig. and ®. 
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FIG. 7. — Relative difference between the output temperature estim ated 
through a one-component fit and the column temperature defined in Eq. ji5\ 
for sources B[ and B2 (panels (a), (b)). The output spectral indices are also 
shown in panels (c), (d), where the dotted lines mark the input value. Panels 
(e) and (f): T-/3 relationship. The dashed line marks the estimated output 
value, while the red circles show the input column temperature as a function 
of the input spectral index. The same parameters estimated through a least 
squares fit are also shown for comparison (blue dots). All the histograms are 
normalized to 1. The light and dark purple contours identify the recovered 
correlation within 1-cr and 2-cr error, respectively. 



This is the most challenging case considered, as the SEDs are 
the combination of two modified black-body, having tempera- 
tures just few degrees apart. It is therefore difficult to identify 
the presence of two components through the comparison be- 
tween Mi and M2, especially when the cold component has 
the same column density as the warm one and results, then, 
to be fainter (cases Bl and B2). For these reasons, the major- 
ity of sources is assigned to Ml in those cases, even if they 
are constituted by two emitting sources. In cases B3 and B4, 
since the cold component is brighter by construction, it is eas- 
ier to detect the presence of two sources, even if the majority 
of sources is still assigned to Mi and many sources have bad 
fit. Where the presence of the second component is more ev- 
ident (cases B3 and B4), the one-component approximation 
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FIG. 8. — Relative difference between the output temperature estim ated 
through a one-component fit and the column temperature defined in Eq. H5i 
for sources B3 and B4 (panels (a), (b)). The output spectral indices are also 
shown in panels (c), (d), where the dotted lines mark the input value. Panels 
(e) and (f): T-/9 relationship. The dashed line marks the estimated output 
value, while the red circles show the input column temperature as a function 
of the input spectral index. The same parameters estimated through a least 
squares fit are also shown for comparison (blue dots). All the histograms are 
normalized to 1. The light and dark purple contours identify the recovered 
correlation within 1-cr and 2-cr error, respectively. 



is clearly worst and the parameters are recovered with larger 
uncertainties. In panels from (a) to (d) of Fig. (|8) we show 
the relative differences between the input and output column 
temperatures and spectral indices for these cases. Cases Bl 
and B2 are affected from the same systematic bias observed 
in the ISM simulations, due to large calibration uncertainties. 
Their parameters are therefore shifted of ~ 15% of the origi- 
nal value, on average. This creates also a spurious T-/3 anti- 
correlation which is nonetheless negligible within 3-ers. The 
same is not true for cases where the cold component is denser 
and the column density therefore higher and comparable with 
the column density of the envelope. Here again the parame- 
ters, mostly the spectral indices, are affected by some bias due 
to the wrong modeling spectral shape. Nonetheless, the over- 
all T—j3 relationship is correct within the uncertainties and is 
particularly accurate for case B4, where the presence of the 
850^m point better constraints the RJ slope and, therefore, 
the spectral index. The T— /3 relationships are presented in 
both figures in panels (e) and (f). 
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An application of the method on a set of cold cores on an 
inter-arm region of the Galactic Plane observed in the Her- 
schel survey, is reported in Sec. (0. 

4.3. HII regions and dust 

Results for HII regions-type of sources are shown in Fig. (|9) 
and reported in Tab. (0. 
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FIG. 9. — Comparison between output and input parameters for sources Ci 
and C2. All the histograms are normalized to 1. Panels from (a) to (d): 
relative difference on the cold and warm temperatures recovered with M2 
with respect to the input values. Panel (e): T-/3 relationship from the sources 
Ci identified as single components. The red circles mark the input column 
density while the light and dark purple contours identify the recovered corre- 
lation within 1-cr and 2-a error, respectively. The same parameters estimated 
through a least squares fit are also shown for comparison (blue dots). 



In both cases the presence of two components is well detected 
with the model identification, especially when the 24/im flux 
is present to constrain the peak of the warm component (case 
C2). As already outlined in Sec. [3] in this analysis we con- 
sider the 24/im flux to be dominated by BG. For this kind of 
sources, M2 fits better than Mj in 100% of cases. When ex- 
cluding the MIPS point (case CI) and we sample the SEDs 
from 70/im down to 500/im, ~ 70% of cases are properly 
classified. Panels from (a) to (d) show the correct recovery 
of the two component temperatures (ISM and HII) by fitting 
the SEDs with M2. They are recovered with a maximum rel- 
ative error of 50% in the worst cases (panels from (a) to (d)). 
Panel (e) shows the T-/3 relationship for sources CI, by fit- 



ting the SEDs from 70/im to 500/im with a modified single 
black-body. The 70/im flux might still be affected from the 
warm component emission and this creates an average shift 
towards warmer temperatures in the parameters recovery. De- 
spite a slight anticorrelation is detected on the data points, it 
results to be negligible when the calibration uncertainties are 
taken into account. 

5. STARLESS CORES IN THE HERSCHEL/HI-GAL SURVEY 

After having tested the validity of the method on simulated 
datasets, we apply it on a s et of starless cores det ected in the 
Herschel Hi-GAL survey dMolinariet al.l 120 1 Qallbb in an m- 
terarm region of the Galactic plane. Hi-GAL is one of the 
Herschel open time key projects to map the entire Galactic 
plane in 2° x 2° tiles with PACS and SPIRE in parallel mode, 
in the bands 70, 160, 250, 350, 500 fj.m. Once the survey 
will be complete, the sky coverage will be 0° < I < 360° and 
-1° < I < 1°. The dataset, due to the spectral range, sensi- 
tivity and sky coverage, is particularly sensitivity to the very 
early stages of high-mass star formation. In this paper, we fo- 
cus on one of the Science Demonstration Phase (SDP) fields 
of the survey: a 2° x 2° tile centered on {I, b) = (59°, 0°). 
This field has been observed in November 2009 and, since 
then, it has been largely studied and used to test the analy- 
sis algorithms. The tile covers an interarm region, since the 
line of sight is tangent to the Sagittarius arm, at an heliocen- 
tric distance between 2 and 7 kpc, where most of the sources 
are located. The overall star formation activity is mostly 
do minated by the Vul pecula OB association (see for exam- 
ple, iBillot et all (120101) '). a molecular complex with triggered 
star forming activity taking place. Nonetheless, the Star For- 
mation Rate (SFR) in that region, estimated from the 70/im 
colors of the observed Young Stellar Objects (YSOs) is low 
(2.6 x 10" 6 M^ /vr. IXfeneziani et al.l dMTl 1. Since the back- 
ground emission is also weak, the source detection is very 
accurate down to very low signal-to-noise thresholds and this 
makes it a favorable set to study molecular cold clumps in a 
prestellar evolutionary stage. 

5.1. Observations and data analysis 

In this context, we us e the data from the Spitzer legacy sur- 
vey MIPSGAL at 24/^m lCarey et al.1 (120091) and the Herschel- 
HiGAL data at 70, 160, 250, 350, 500 /im. Kine matic sun dis- 
tances are also estimated from CO observations iRusseil et alj 
(12011 . 

Herschel maps have been produced by means of the RO- 
MAGAL algorithm, based on a Generalized Least Square 
technique, after a careful preprocessing in which glitches and 
systematic effects present in the data have been removed. For 
further information about the pipeline follow ed from raw data 
up to s ky maps delivery, we refer the reader to lTraficante et ID 
(120111) . The source detection and extraction has been per- 
forme d by means of the CuTEX algorithm iMolinari et al.l 
(1201 ll) which double-differentiates the sky image and study 
the curvature variations above a given threshold. The identi- 
fied source profiles are then fitted with a 2D Gaussian plus an 
underlying inclined planar plateau. The Herschel catalog is 
then produced in the following way: when a source is found 
in the longest wavelength map (500/im), we look for an asso- 
ciation in the following band (350/im) in a radius as large as 
the FWHM of the largest beam between the two. If a source is 
found, then we proceed to the association with the following 
band (250/im) and so on up to the 70/im map. If more than a 
source is found in the radius, we associate the nearest. Fluxes 
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corrections are applied both in the PACS and SPIRE bands in 
order to take into account of the non perfectly Gaussian shape 
of the beam. By merging the Herschel catalog with the MIPS 
catalog by means of the same procedure described above, we 
have an overall catalog of the entire region from 24 to 500 
fjm. 

The procedure to select and identify the starless molecular 
clumps is summarized in the following: 

1. the sources are not detected either in the 70/im nor in 
the 24/im band in order to exclude the presence of a 
forming star inside the clump; 

2. the sources are detected in all the four remaining bands, 
from 160/iin to 500/im. This makes us sure not to in- 
clude spurious detections and to have a better charac- 
terization of the SEDs; 

3. the sources are gravi tationally bound. As 
in IVeneziani et al.l ( 12013b . we identify as bounded 
objects the sources with mass M ^ 0.5 Mbe where 
Mbe is the Bonnor-Ebert mass. A Bonnor-Ebert sphere 
is an isothermal sphere at hydrostatic equilibrium. 
Therefore, in absence of internal turbulence and 
assuming thermal pressure, the BE mass is a good 
approximation of the virial mass. 

The absence of turbulence is a necessary approximation due 
to the fact that no spectroscopic data are available on this sam- 
ple of sources. We cannot then check the internal motions 
occurring in the clumps and see if there are supporting mech- 
anisms other than thermal pressure. After all these criteria 
have been applied we have a total of 103 sources. Since only 
four fluxes are available, we assign an upper limit of 0.2 Jy 
to the 70/iin band, in order to perform the model identifica- 
tion. With Hi-GAL sensitivity and coverage, 0.2 Jy is the 
minimum flu x detected in the PAC S blue band even at a 15 
kpc distance (IVeneziani et al.l 120131) . 5 sources are excluded 
through the model identification, leaving us with a sample of 
91 sources. The majority of those objects have radius 5 > 0.1 
pc, confirming that they are essentially clumps. 

5.2. SEDs fitting and dust parameter estimate 

We associate a conservative error bar of 30% of the fluxes 
to the SEDs from 160/im to 500/iin. This error takes into ac- 
count the gaussian random noise, the source multiplicity, the 
background fluctuations and the overall uncertainty on fluxes 
recovery of the extraction and photometry algorithms. The 
calibration errors included are PACS and SPIRE uncertainties 
on extended emission, the same we considered in our simu- 
lation, i.e. 20% for the 160/^m PACS band and 15% for the 
250, 350 and 500^ SPIRE bands. 

The Bayesian fit with multi-variate Gaussian priors and fol- 
lowing marginalization on the calibration uncertainties are 
then run to estimate the SEDs parameters. After a first MCMC 
run with uniform wide priors, we fix the center of the multi- 
variate Gaussian priors in the point [To, log(eo), A)] = [12, 
-9, 2]. As already described, even if the full sample is col- 
lected in the same observation, the covariances are estimated 
source by source since they depend on random fluctuations, 
SED sampling and overlap of other sources along the line of 
sight, assuming that all the other sources of error, i.e. sys- 
tematic effects, are the same for all the objects in the same 
observation. In each covariance matrix, the standard devia- 
tions are enlarged by ten times in order to make sure not to 
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FIG. 10. — Example of SEDs of four cold clumps at different temperatures 
and spectral indices values. The black symbols mark the fluxes with the as- 
sociated 30% error bars. The blue dashes indicate the calibration uncertainty. 



constrain the final parameters with a too tight distribution but 
keeping the same correlations. The values are nevertheless 
well described by the ones obtained from the simulation run 
of cold cores and reported in Tab. (3). 

An example of the SEDs with calibration and statistical un- 
certainties and their best fits, is shown in Fig.flOl 

The final temperature and spectral index distributions, as 
well as the T-/3 relationship recovered with our bayesian 
method (BM), are shown in Fig. ( fTTT i as black dots with the 
associated error bars. Their median values are (^59) = 11.8 
K and (^59) = 1 .6. The relationship between the temperature 
and the spectral index has parameter values Afjg = 1 .9 1 ± 0.29 
and a™ = 0.28 ± 0.29. In order to compare the results of 
our Bayesian method with a technique based on least squares 
minimization, we fit the SEDs by associating to the fluxes 
only a statistical 30% error bar, without taking into account 
the calibration uncertainties. The temperature and spectral in- 
dices values are shown in Fig. ( fTTT i as red dots. Results from 
the least square fitting show a spurious clear anticorrelation: 
A" 9 = 0.8 ±0.3 and a" 9 = -1.3 ± 0.5 which is, nonetheless, 
consistent with no anticorrelation within a 3er range. 

As already described in the Introduction an increase of 
the emissivity spectral index in cold environments had been 
observed both in la b experimen t s dMenv et ail 120071) and 
in real observation (iDupac et al.l [2003b lDesert et al.l 120081 
IVeneziani et alJ2010tlParadis et al.l20ldlBracco et al.ll201 lT) . 
It can be explained in two ways: t he spectral index is tempera- 
ture dependent, as suggested from lMenv et al. (2007) from lab 
experiments on amorphous silicate-based dust grains, or the 
spectral index increases as a consequence of grain aggregation 
in col d and dense environments (see for example [Ormel et aD 
1201 11) . In the observed Herschel cold cores, we do not detect 
a significant increase of the spectral index with respect to the 
ISM, since j3 values are spread between 1 and 2.5 (Fig. QT| 
middle panel) and, when the T-/3 trend is studied, points are 
scattered on the T-/3 plane excluding the possibility of a corre- 
lation (Fig.QT| right panel). The absence of a correlation and, 
more in general, the lack of an increase of the spectral index 
with respect to the diffuse ISM, can be explained considering 
that, as pointed out several times during the paper, this anal- 
ysis is very sensitive to systematic effects which might range 
from instrumental issues other than calibration, (i.e. photo- 
metric extraction, catalog band-merging), to multiple over- 
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lap of sources along the line of sight. A possible physical 
explanation of our findings might be that we are observing 
cores of differen t ages and this reflects the spre ad of the spec- 
tral indices. As lOssenkopf & Henningl (119941) shows, when 
silicate dust grains are covered by a thick ice mantle, their 
emission becomes independent of the optical properties of the 
material underlying the mantle, as to mask the T-/3 anticor- 
relation. Therefore, depending on the mantle thickness, i.e. 
the timescale of the core, the optical properties of the grains 
and th eir spectral indices change. Moreover, as lOrmel et"aD 
(120091) shows, depending on the cold core timescale, grains 
might or might not have the time to aggregate and this also 
can explain the spread of the spectral index values. 

6. CONCLUSIONS 

We have presented a method, based on Bayesian statistics, 
to fit for intrinsically correlated parameters and look for the 
underlying correlation law. This method is different from the 
ones previously used for similar analysis, as we make use of 
Bayesian statistics, which fully sample the parameter space 
taking into account the covariance among different parame- 
ters, together with a proper treatment of systematic errors by 
means of a Monte Carlo procedure. We apply the method 
to estimate the physical properties and the relationship be- 
tween the temperature and the spectral index of a set of simu- 
lated astrophysical sources detectable in the Herschel PACS 
and SPIRE bands (between 70/im and 500/um), including 
also Planck-HFI (850^m), IRAS (100 /im) and MIPS (24/im) 
bands, in order to better constrain the colder and warmer com- 
ponents. The sources are chosen over a wide range of temper- 
atures and compositions in order to test the efficiency of the 
method with different samplings of the SEDs. For this pur- 
pose, we consider the cases of a one-component ISM (with T f / 
and j3 either correlated, case A2 and A3, or uncorrelated, case 
Al), and of two-component temperature sources both warm 
and cold. The case of two-comp onent sources is particu larly 
interesting because, as shown in IShetty et"aT1 (l2009albl) . the 
combination of more than one source along the line of sight 
generates a spurious anticorrelation by itself. 

The most relevant results are: 

• The input physical values and their correlations are 
well recovered in the Herschel bands when we observe 
an ISM like sample both without and with correlation 
(cases Al, A2 and A3, respectively); 

• when two sources just few degrees apart are combined, 
as cold clumps, it is difficult to identify the underly- 
ing model especially when the core and the envelope 
have the same column density (cases Bl and B2) and 
the core emission is therefore much fainter. If the core 
is thicker it is easier to detect (cases B3 and B4) and, in 
this case, the presence of the 850/im flux (B4) helps to 
better constrain the RJ part of the SEDs and, therefore, 
the physical parameters and their correlations. Even if a 
spurious anticorrelation is detected in the cold clumps, 
when taking into account the calibration uncertainties, 
the detected anticorrelation is negligible within 3a. 

• when the source is a combination of two components 
several degrees apart, as for HII regions, the underly- 
ing model is correctly identified, especially when the 
24/im flux is present. Also in this case, a slight anti- 
correlation is measured in the dust component but the 



detection is neglibile when taking into account calibra- 
tion uncertainties. 

The method has been applied to a sample of starless clumps 
in a 2° x 2° field centered on (£, b) = (59°, 0) observed in 
the Herschel/Hi-GAL survey. Being an inter-arm region of 
the Galactic plane, it is not not very dense with sources and 
the star formation activity is therefore very low. The average 
temperature and spectral index of the sample are ~ 11.8 K 
and 1.6, respectively, as expected for cold sources. We do not 
detect a T-/3 correlation within the uncertainties. 
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FIG. 1 1 . — Physical parameters of the starless clumps in the Herschel/Hi-GAL (£, b) = (59°, 0) field. First and second panels show the temperature and spectral 
index distributions, respectively. The third panel shows the T-/3 distribution. The light and dark purple contours identify the recovered correlation within l-cr and 
2-cr error, respectively. The red dots mark the temperatures and spectral indices obtained by fitting the SEDs with a least square method, for comparison. 
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